#include <iostream>
#include "includes/dg-evolve-finite-element.hh"
// #include "includes/basis-function.hh"
DGEvolveFiniteElement::DGEvolveFiniteElement(){
};
DGEvolveFiniteElement::~DGEvolveFiniteElement(){
};
void DGEvolveFiniteElement::dg_evolve_finite_element(SubprojectBase* subproject,double* globalcoord, double* unitouternormal,double pillarlen,int i,int j, int currentdim, int* interscounter, double rdt){
  double var_i[memory::phynumber];
  double fequ_i[memory::phynumber];
  double dbasis[memory::dgorder];

  double var_i_left[memory::phynumber];
  double var_i_right[memory::phynumber];
  double var_i_middle[memory::phynumber];      
  double fequ_i_left[memory::phynumber];
  double fequ_i_right[memory::phynumber];
  double fequ_i_middle[memory::phynumber];
  double dbasis_left[memory::dgorder];
  double dbasis_right[memory::dgorder];
  double dbasis_middle[memory::dgorder];  

  basis.getVarLeft(i,var_i_left);
  basis.getVarRight(i,var_i_right);
  basis.getVarMiddle(i,var_i_middle);
  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_i_left,fequ_i_left);
  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_i_right,fequ_i_right);
  subproject->equation->get_equation_flux(gridsetting::dimension,unitouternormal,var_i_middle,fequ_i_middle);
  basis.getDBasisLeft(dbasis_left);
  basis.getDBasisRight(dbasis_right);
  basis.getDBasisMiddle(dbasis_middle);
  
  // // trapezoid quadratureu
  // for(int tmpj=0; tmpj<memory::dgorder; tmpj++)
  //   for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
  //     memory::update[timerule::roundnumber][tmpi*memory::phynumber+tmpj][i]=pillarlen*(fequ_i_left[tmpi]*dbasis_left[tmpj]/2.0+fequ_i_right[tmpi]*dbasis_right[tmpj]/2.0);

  // simpson quadrature
  for(int tmpj=0; tmpj<memory::dgorder; tmpj++)
    for(int tmpi=0; tmpi<memory::phynumber; tmpi++)
      memory::update[timerule::roundnumber][tmpj*memory::phynumber+tmpi][i]=pillarlen*(fequ_i_left[tmpi]*dbasis_left[tmpj]/3.0+fequ_i_middle[tmpi]*dbasis_middle[tmpj]*4.0/3.0+fequ_i_right[tmpi]*dbasis_right[tmpj]*3.0);

  if(i==0)
    std::cout << "dg evolve for finite element" << std::endl;

};
